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1 INTRODUCTION 



Gamma-ray burst (GRB) afterglows can be explained from 
the interaction between an initially relativistic shock wave 
of hot fluid and the medium surrounding the burster. On 
passage of the shock electrons get accelerated to relativis- 
tic velocities (even with respect to the already relativistic 
local fluid flow) and small scale magnetic fields are gener- 
ated. Under influence of the magnetic field, the electrons will 
produce synchrotron radiation, which will be seen by the ob- 
server. This model has been very successful when applied to 
broadband afterglow data, but thus far model predictions 
have been made using simplifying assumptions for the blast 
wave structure (approx imating the blast wave width by a 
homogeneous slab, e.g . Wiiers et al.l 19971 ; iMeszaros fe Reesl 



ABSTRACT 

We present a study of the intermediate regime between ultra-relativistic and non- 
relativistic flow for gamma-ray burst afterglows. The hydrodynamics of spherically 
symmetric blast waves is numerically calculated using the AMRVAC adaptive mesh re- 
finement code. Spectra and light curves are calculated using a separate radiation code 
that, for the first time, links a parametrisation of the microphysics of shock accelera- 
tion, synchrotron self- absorption and electron cooling to a high-performance hydrody- 
namics simulation. For the dynamics we find that the transition to the nonrelativistic 
regime generally occurs later than expected, that the Sedov- Taylor solution overpre- 
dicts the late time blast wave radius and that the analytical formula for the blast 
wave velocity from iHuang et al.l (|1999| ) overpredicts the late time velocity by a factor 
4/3. Also we find that the lab frame density directly behind the shock front divided 
by the fluid Lorentz factor squared remains very close to four times the unshocked 
density, while the effective adiabatic index of the shock changes from relativistic to 
nonrelativistic. For the radiation we find that the flux may differ up to an order of 
magnitude depending on the equation of state that is used for the fluid and that the 
counterjet leads to a clear rebrightening at late times for hard-edged jets. Simulating 
GRB030329 using predictions for its physical parameters from the literature leads to 
spectra and light curves that may differ significantly from the actual data, emphasizing 
the need for very accurate modelling. Predicted light curves at low radio frequencies 
for a hard-edged jet model of GRB030329 with opening angle 22 degrees show typi- 
cally two distinct peaks, due to the combined effect of jet break, non relativistic break 
and counterjet. Spatially resolved afterglow images show a ring-like structure. 

suiting spectra (see lGranot et al"1l200ll ; Ipownes et al.| [2002 V 
More recent simulations have been used to address the spe- 
cific theoretical issue of the visible effect of the blast wave 



homo geneous slab, e.g. Wuer s et al.1199 /; Mcsz aros fe Kee; 
Il997l ; ISari et al.lll998l ; TRhoadsj|l999l ) or from analytical so 



lutions in e ither the ultrarelat i vistic or the nonrelativistic 
regime (e.g. Granot et al.ll 19991; iGruzinov fc Waxmanll 19991 ; 
IWfiers fe Galamalll999l : lFrail et alJl200Ch . 

Since the beginning of the decade, fluid simulations have 
been performed to study afterglow blast waves and their re- 



encountering a density perturbation (Nakar & Granot 2007: 
Ivan Eerten et al.ll2009h . Very recentlv lZhang fe MacFadvenl 
12009) studied the transition to the transrelativistic regime 
and the spreading of a collimated outflow, using an adaptive 
mesh technique for the fluid simulation. They made some 
simplifying assumptions for the radiation mec hanism, when 
comp ared to the early analytical efforts (e.g. iGranot et al.l 
1 1999h . such as approximating the cooling time by the lab 
frame time and ignoring synchrotron self-absorption. 

The aim of this paper is to present a theoretical and 
qualitative study of the transition regime between relativis- 
tic and nonrelativistic blast waves and the effect on the light 
curves and spectra at various wavelengths, using adaptive 
mesh relativistic fluid simulations for blast waves from an 
explosion in a homogeneous medium, while including all 
details of the synchrotron radiation mechanism that have 
been used for earlier analytical estimates. Also we present 
resolved afterglow images. We study spherical blast waves 
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and sharp edged jets obtained by taking conic sections from 
a spherically symmetric fluid flow. 

Obviously these simulations do not yet fully address the 
complete GRB afterglow picture of a realistic, 2D-dynamical 
jet, which we address in future work. However, some GRB 
afterglows have power law decays that last for months with- 
out a jet break, and thus may be (nearly) spherical. These 
are of course already addressed in the present work. Also, by 
studying the conic sections from spherical flows, we already 
address some aspects of jet behaviour, which allows us to 
probe some outstanding issues, such as whether the reced- 
ing jet may lead to visible features in the late light curve, and 
whether a dynamical jet break must be truly achromatic. Fi- 
nally, any fluid flow behaviour typical to higher-dimensional 
simulations, like lateral spreading of the jet, is best under- 
stood from a direct comparison to one-dimensional simula- 
tions and its effects on the light curve will in practice be 
modeled as a deviation from the heuristic description based 
on analytical approximations and one-dimensional simula- 
tions (i.e. as an additional smooth jet break). A compan- 
ion paper is in preparation that will discuss the practical 
consequences for broadband afterglow data fitting from the 
underlying model from this paper. 

This paper is organized as follows. In section [2] we dis- 
cuss our radiation code and how it expands upon an ap- 
proach outlined earlier in I Van Eerten fc Wiiersl (|2009l ), here- 
after EW09. A proper treatment of synchrotron radiation 
and shock wave generation of accelerated particles and small 
scale magnetic fields requires us to trace some additional 
quantities along with the fluid quantities. 

In section [3] we provide the details of our simulations 
that assume typical GRB parameters. We show how the 
blast wave starts out in the ultrarelativistic regime and 
smoothly approaches the nonrelativistic regime. We discuss 
the consequences of different equations of state for the fluid 
and how our simulations differ from analytical approxima- 
tions for the nonrelativistic regime. We show how the fluid 
lab frame density divided by the fluid Lorentz factor squared 
right behind the shock remains always close to four times 
that in front of the shock, even though we have differing 
adiabatic indices in both the relativistic and nonrelativis- 
tic regimes. Three additional quantities needed to be traced 
and we present results for the behaviour of these three: the 
accelerated electron number density, the magnetic field en- 
ergy density and the accelerated particle distribution upper 
cut-off Lorentz factor. We explain how calculation of the lat- 
ter especially is numerically challenging and how it shapes 
the spectrum beyond the cooling break. 

In section U we take our results from section [3] and cal- 
culate spectra and light curves. We calculate spectra at 1, 
10, 100, 1,000 and 10,000 days in observer time. We sepa- 
rately discuss the different factors contributing to the shape 
of the light curves: the equation of state, the evolution of the 
magnetic field and the evolution of the accelerated particle 
distribution. 

We then turn to the specific case of GRB 0303029 in sec- 
tion [S] We take the explosion parameters that have been es- 
tablished for this burst by previous authors to set up a simu- 
lation. We qualitatively compare the resulting light curves to 
radio data at different wavelengths, assuming both a spher- 
ical explosion and a hard edged jet with opening angle of 
22 degrees. We provide spatially resolved radio images and 



make a qualitative prediction for the expected signal at radio 
wavelengths that will be observable with the next generation 
of telescopes, like lofar. 

We discuss our results in[S] In the appendices we pro- 
vide additional technical details on the numerical implemen- 
tation of our approach and a discussion on the theoretical 
limitations and assumptions of our approach. 



2 THE RADIATION CODE 

In this paper we follow the approach first outlined in EW09, 
where we calculate spectra and light curves from the output 
of a relativistic hydrodynamics (RHD) code using a separate 
radiation code. For the RHD simulations we use AMRVAC, a 
high performance code that includes adaptive-mesh refine - 
ment (AMR) fsee lKeppens et al.|[2003l : iMeliani et al.ll2007r i. 
AMRVAC calculates the evolution of the following conserved 
variables: 

D = "/p', S = r y 2 h'v, t = -y 2 h' — p' — -yp' c 2 , (1) 

with 7 the Lorentz factor, p' the proper density, h! the rel- 
ativistic (i.e. including rest mass) enthalpy density, v the 
three velocity, p' the pressure and c the speed of light. In 
the entire paper, all comoving quantities will be primed. 

In the second stage we use a radiation code to obtain 
the received flux for a given observer frequency, time and 
distance, from the local values of conserved variables at any 
contributing point in the fluid (we also use two auxiliary 
quantities, j and p' , that AMRVAC stores as well in order 
to facilitate its calculation of the time evolution of the con- 
served variables). The radiation mechanism that is consid- 
ered is synchrotron radiation and a number of parameters 
have been introduced in EW09 that capture the underlying 
radiation and shock microphysics. There are four of these 
'ignorance' parameters. The fraction of the thermal energy 
that resides in the tangled-up magnetic field that is gen- 
erated by the passage of a shock eb usually has a value 
around 0.01. The fraction of electrons £n that is accelerated 
into a relativistic power law distribution in energy also by 
the passage of a shock is usually of order unity in the rel- 
ativistic regime. The thermal energy fraction captured by 
these electrons ee ^ 0.1 and minus the slope of the electron 
distribution p ^ 2.5. 

The flux calculated by the radiation code is given by 

F v = l -±^ [ -^(1 - McdA&U. (2) 

Here z denotes redshift, r Q b s denotes the observer luminosity 
distance, d 2 Pv/di^df2 the received power per unit volume, 
frequency and solid angle, dA the equidistant surface ele- 
ment given by the intersection of the fluid grid with that 
surface from which radiation is poised to arrive exactly at 
fobs, and t a the emission time. Note that in this terminology, 
flux is defined per unit frequency. The integral 

J (1 - Pn)c dAdt c 

is effectively an integral over the entire radiating volume. 
p is the angle between the local fluid velocity and the ob- 
server position, P the fluid velocity in units of c and the 
factor (1 — /3p) is a retardation effect due to the moving 
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of the radiating source. The detailed dependency of the re- 
ceived power on the ignorance parameters and local fluid 
conditions is explained in EW09. However, in that paper 
only ultra-relativistic flows were addressed and in order to 
include subrelativistic and nonrelativistic flows as well, a 
number of features were added to our radiation code. Also 
we have added synchrotron self-absorption and the possibil- 
ity to resolve the signal from the fluid into an image on the 
sky. We now have a generic radiation code that is capable 
of calculating the spatially resolved synchrotron radiation 
profile from an arbitrary fluid flow. The additional physics 
that we have included is explained below, with some of the 
practical numerical issues discussed separately in appendix 

m 



2.1 Realistic equation of state 

In EW09 we applied a fixed adiabatic index F a d equation of 
state (EOS) 



o = (r ad - l)e' th , 



(3) 



where e' th is the thermal energy density. In practice r a d was 
always set to 4/3. However, when following a fluid from the 
relativistic regime (with flow velocities ^ c and thermal en- 
ergy density dominating the rest mass energy density) down 
to the classical regime, this fixed adiabatic index b ecomes 
too r estrictive. We therefore apply a Synge-like EOS ijSvngd 
Il957f ) that results in an effective adiabatic index varying 
smoothly from 4/3 to its classical limit 5/3: 



P'c 2 



P'c 2 



p'c 2 



(4) 



where e' denotes the comoving energy density including rest 
mass, e' = p'c 2 + e[ h . This EOS has al ready been applied 
in amrvac Csee lMeliani et al.l f2008. 200|). Also, because the 
radiation code reads both the conserved variables as well as 
p' from disc directly, it does not invoke any EOS itself, and 
no change in the radiation code was needed. The resulting 
effective adiabatic index is given by 



r 5 1 (^ P' 2(A 

1 ad, eft = — — — 1 ST" 



(5) 



The effect of an advanced EOS on the behaviour of the fluid 
is profound and we discuss this in detail in section 14.31 



2.2 Electron cooling 

The shape of the observed spectrum from a single fluid 
cell, if electron cooling does not play a role, follows directly 
from the dimensionless function Q(y' /v m ), first introduced 
in EW09. It has the limiting behaviour Q oc (y' /u^) 1 ' 3 for 
small (i//i4) and Q oc (^'/^) (1 " p)/2 for large (u 1 ' /v' m ). The 
received power depends on this shape and on the local fluid 
quantities via 



AvAQ. 



(p - l)V3g* &nB' 
8irm e c 2 7 3 (1 — Pp) 3 



7> ■ (6) 



Here n denotes the lab frame number density (of all elec- 
trons, both accelerated and thermal). B' denotes the local 
comoving magnetic field strength, calculated from the ther- 
mal energy density after the passage of a shock. m c and 



q e denote electron mass and charge, respectively. The fre- 
quency v' m is the synchrotron peak frequency, and it is re- 
lated to the lower cut-off Lorentz factor 7 m of the power law 
accelerated electrons via 

/ 3<? c 



(7) 



47rm e c 

If cooling plays no role, the evolution of 7 m is completely 
adiabatic, which has as a consequence that the total fraction 
of the local thermal energy density residing in the power-law 
accelerated particle distribution remains fixed. y' m will be 
related to e' th throughout the downstream fluid according to 



7m = 



^n'm e c 2 



(8) 



When cooling does play a role, however, this picture is 
changed. It now becomes necessary to introduce an upper 
cut-off Lorentz factor 7 M as well. In a single fluid element, no 
accelerated electrons with energies above 7 M will be found, 
because these have cooled to energies at or below 7 M . The 
temporal evolution of any electron Lorentz factor 7^, and 
therefore of 7 m and 7m as well, is given by 



d 7c 
At' 



3n' 



An' 



- «) 



■ a T B r 



(9) 



At' 6nm e c 
where ctt is the Thomson cross section, m e the electron mass 
and t' the comoving time. The final term in this equation 
reflects synchrotron radiation losses, and if it is omitted only 
the adiabatic cooling term is left and it can be shown that 
this will result in the aforementioned fixed relation between 
7e and e' th . In light of the previous subsection on the EOS, 
it may be worth noting that equation [5] is derived for a 
relativistic electron distribution with adiabatic index 4/3 
and that this remains valid even if the bulk of the fluid 
becomes nonrelativistic. After all, the power-law accelerated 
electrons are relativistic by definition. 

The above has the following consequences for the sim- 
ulations and radiation code. Because for low values of 7e, 
the radiation loss term can be neglected next to the adia- 
batic expansion term, we will not apply equation [§] to 7 m 
and we will continue to calculate 7 m locally using equation 
[8] in the radiation code. For 7 M this is not an option and we 
numerically solve equation in AMRVAC, resetting 7 M to a 
high value wherever we detect the passage of a shock. This 
reset implements the shock-acceleration of particles. In ap- 
pendix [A] we discuss the numerical issues of this approach 
in some more detail. Also, in appendix IB2I we show that the 
gyral radius even for the high energy electrons contributing 
to the observed spectrum (within the frequency range under 
consideration, 10 8 — 10 18 Hz) is orders of magnitude smaller 
then the relevant fluid scales. 

Finally, we summarize the consequences of electron 
cooling for the spectrum discussed earlier in EW09. The 
received power from a fluid element is now given by 



d 2 P v 



oc 



£vinB' 



-Q 



(10) 



AvASl " 7 3 (1 - finf ' 

with the relations between the upper and lower cut-off 
Lorentz factors and their corresponding critical frequencies 
given by equation [7] Q is a generalisation of Q and the flux 
at frequencies v' above u' m drops exponentially. That the 
resulting spectrum from the entire fluid does not show an 
exponential drop is due to the fact that there will always be 
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some fluid elements contributing for which v' m is still suffi- 
ciently high. The effect of this 'hot region' close to the shock 
front (with a size that depends on the observer frequency) 
on the composite synchrotron spectrum from a shock will 
be a steepening of the slope by —1/2 instead. The cooling 
break is found at that frequency for which the width of the 
hot region becomes comparable to the width of the blast 
wave. 



2.3 Magnetic field energy evolution 

The magnetic field directly behind a shock has been 
parametrised using 



ese th . 



(11) 



Furthermore we assumed the number of magnetic flux lines 
threading a surface comoving with a fluid element to remain 
invariant, resulting in 



e B oc P 



14/3 



(12) 



For relativistic fluids this implies that the fraction es re- 
mains fixed downstream, because e' th oc p' rad . For a chang- 
ing adiabatic index it is no longer possible to calculate e B 
a posteriori from e' tb , since the relation between the two is 
now no longer fixed. It becomes necessary to numerically 
solve in AMRVAC the equation 



At' p' 4 / 3 



0. 



(13) 



Like 7m, we reset e' B whenever a shock is encountered. The 
practical implementation of the evolving magnetic field is 
again discussed in appendix [S] We note here that the as- 
sumption of frozen field lines is not essential, and that we 
can in principle include different magnetic field behaviour 
either by adding a source term to equation [T3] (parametris- 
ing for example, magnetic field decay through reconnection) 
or by implementing a different equation entirely. 



2.4 Changing fraction of accelerated particles 

Although £jv, the fraction of electrons accelerated by the 
passage of a shock is often assumed to be of the order unity 
for highly relativistic blast waves, it has to be lower at late 
times because otherwise there would not be enough energy 
available per accelerated electron to create a relativistic dis- 
tribution (in other words, to ensure that 7^ > 1). We have 
implemented this change in our code by replacing user pa- 
rameter £at by £n,nr, that is, the fraction of electrons that 
is accelerated in the nonrelativistic limit. The fraction at 
the relativistic limit we set to one. Because 7/? is the most 
direct measure of how relativistic the fluid flow locally is, we 
have parametrised the simplest possible smooth transition 
between both limiting cases by 



f3j_ + £n,nr. 
1.0 + 07 ' 



(14) 



Whenever the passage of a shock is detected, AMRVAC resets 
the number density of accelerated electrons n^ c according to 
"■ace = £n*i', with £jv determined using the equation above. 
As with the magnetic field energy density, we now need to 



follow n' acc explicitly. Because n' acc is a number density, its 
evolution is described by a continuity equation, following 

^"-acc7 + ■^-n'^civ 1 = 0, (15) 
and is therefore easily implemented in AMRVAC. 



2.5 Synchrotron self-absorption 

In previous work we have solved equation [2] by first integrat- 
ing over A for a given emission time t e (and thus for a single 
snapshot), followed by an integration over t c . If we switch 
the order of the integrations then the integral over t e rep- 
resents the solution to a linear radiative transfer equation 
without absorption, with the intensity given by 



AuAQ 



(l-^)cdte 



(16) 



The integral over A then represents a summation over all 
rays. The full linear radiative transfer equation including 
synchrotron self-absorption has the form 



AI V 

— — = -a v l v + >, 
Az 



(17) 



with j v = A 2 Pv/ AvAQ, and dz = cAt b s = (1 — f3fj,)cAt e 
The synchrotron self-absorption coefficient is given by 



7m 



8-Km c v' 



AP' <e> ,2 

7c 



Av' 



7m 



d%V 7e 2 



d 7c . (18) 



Here AP' <e> / Av' denotes the emitted power per ensemble 
electron and A r g(7e) the electron number density for rela- 
tivistic electrons accelerated to y' e . Integrating A r e(7e) over 
possible electron Lorentz factors yields n' &cc by definition. 
These quantities are defined and explained in detail in EW09 
(see also appendix I A3|) . 

In this treatment of the self-absorption coefficient we 
only take into account transitions between already occupied 
energy levels of electrons, leading to the integration limits of 
equation [18] being exactly 7^ and 7 M . In this way we ignore 
stimulated emission arising from a population inversion be- 
low ~/' m . This results in values of the absorption coefficient 
that are larger by a f actor of 3(p + 2)/4 when compared to 
iGranot fc Sari l|2002h . 

In our radiation code, we now calculate the linear ra- 
diative transfer equation for each individual ray by not inte- 
grating over the two-dimensional surface A (to get a single 
flux value from the collection of rays) until after the final 
snapshot has been processed. In addition to allowing us to 
include the effect of self-absorption, we now also get a spa- 
tially resolved signal from the fluid, showi ng the expected 
ring structure (extending predictions from IGranot fc Loebl 
l200ll to the nonrelativistic regime). We use an adaptive- 
mesh type approach to A in order to ensure an adequate 
spatial resolution, see appendix I A4I 



3 FLUID DYNAMICS 

In this section we describe the setup of our relativistic fluid 
simulations and compare the results against the theoretically 
expected behaviour. 
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3.1 Expected early and late time behaviour 

Both the early and late time behaviour of the fluid can be 
described by a self-similar solution that is determined com- 
pletely from the explosion energy E and the circumburst 
number density no. 

At early stages, the B landford-McKee (BM) solution 
l|Blandford fc McKed[l976T ) for relativistic blast waves pre- 
dicts the following relation between the shock front fluid 
Lorentz factor F and the explosion time (i, which is the 
same as the emission time t e ): 



167rp t 3 c 5 ' y J 

The density po is related to the number density through the 
proton mass: po = m p no. The shock radius R(t) is then 
given by 

ii W = rf ( 1 -l^)- ( 2 °) 

To lowest order R(t) is just ct, while the shock front fluid 
velocity j3 ^ 1. Further analytical equations for the fluid 
profile (in terms of pressure p, Lorentz factor 7, number 
density n, etc.), behind the shock front can be found in 
iBlandford fc McKedll976l . 

At late stages the evolution of the blast wave is de- 
scribed by t he Sedov-Taylor (ST) solution ijSedovl 1 19591 : 
lTavlor|[l950l ). For a fixed adiabatic index 5/3, the shock 
radius is now given by 

/ p.2 \ 1/5 

R(t) « 1.15 (^-j , (21) 

which follows directly from dimensional analysis (except for 
the numerical constant). In this classical approximation, the 
speed of light c does not appear. The shock front Lorentz 
factor is approximately one, while /3 can be found from /3 = 
AR(t) / c&t. Again ana lytical form ulae for the fluid profile 
exist in the literature (|Sedovlll959l ). 

At some point in time the evolution of the blast wave 
will no longer be adequately described by the BM solution 
but will become more and more dictated by the ST solu- 
tion. An estimate for the turning point |Piranl [20051 ) can 
be made by equating the explosion energy to the total rest 
mass energy that is swept up: 

E = p c 2 ^irR^ n . (22) 

Solving for i?NR returns the approximate radius at which 
the original explosion energy in the blast wave is no longer 
dominant over its rest mass energy. 

Analytical estimates for the bulk fluid flow velocity in- 
cluding the in termediate regime a lso exist. One such exam- 
ple is found in lHuang et alJ l| 19991 ) and we discuss it in more 
detail as well as compare their prediction for 7(t) right be- 
hind the shock front directly against our simulation results 
in section [3.3.1l 

3.2 Setup of simulations 

We have performed a number of simulations using the typi- 
cal values for a GRB exploding into a homogeneous medium. 
We set up our simulations starting from the BM solution. 
The isotropic explosion energy E = 1 ■ 10 52 erg, the medium 



number density no = 1cm -3 . We have set the initial shock 
Lorentz factor to 10 (and the fluid Lorentz factor therefore 
7, differing by a factor y/2). Although both AMRVAC and 
the radiation code are able to deal with far higher Lorentz 
factors, the focus for this research is on the transition to 
the nonrelativistic regime and for that purpose this rela- 
tively low Lorentz factor is sufficient. We have continued the 
simulations until the fluid proper velocity in the lab frame 
P - 0.01. 

We have used both the advanced equation of state and 
a fixed adiabatic index at 4/3 and 5/3. In the advanced EOS 
simulation we have also calculated the other quantities men- 
tioned in the previous section: £b,£n and 7^- The value at 
the shock front for eb was set to the standard 0.01 and the 
non-relativistic limit for £n was set at £n,nr = 0.1. Suffi- 
ciently high values for at the shock front are chosen, 
generally on the order of 10 7 . 

In AMRVAC it is the number of refinement levels that 
determines the accuracy of the simulation. We have used 
17 levels of refinement and 120 cells at the lowest refine- 
ment level. The grid was initially taken to run from 10 16 cm 
to 10 19 cm. The effective spatial resolution due to adaptive 
mesh refinement was therefore ^ 1.27- 10 12 cm. This should 
be compared against the width of the blast wave at the start 
of the simulation, when it is the smallest. This is approxi- 
mately equal to R(t) /F 2 ^ 3 • 10 15 cm, for a starting shock 
Lorentz factor of 10. 

Convergence of our results has been checked by per- 
forming simulations at different refinement levels and by 
simulations running for a shorter time on a smaller grid 
(thereby increasing the resolution). For the light curves and 
spectra we have used simulations with a shorter running 
time of 12.2 • 10 3 days. At this stage the fluid velocity di- 
rectly behind the shock is still six percent of light speed, but 
we have full coverage up to 10, 000 days in observer time. 
The corresponding grid size is 1 • 10 17 cm to 6.7 • 10 18 cm, 
leading to an effective resolution of 8.3- 10 11 cm. On a stan- 
dard desktop PC0 the RHD simulations typically took a few 
days to complete and the radiation calculation a few hours. 



3.3 Results 

3.3.1 Blast wave velocity 

The solid line in figure [1] shows ^7 at the shock front for the 
advanced EOS simulation. The expected scaling behaviour 
at the early stage is dictated by F tx t~ 3 ^ 2 and at the late 
stage by /3 oc t -3 ^ 5 . We have plotted this asymptotic be- 
haviour as well, setting the early stage scaling coefficient 
from the initial value at /?7 ^ 7 and the late stage scaling 
coefficient at /?7 ^ 0.016 (this point lies far to the right 
outside the plot). The shock velocity is shown to smoothly 
evolve from the BM solution to the ST solution. The meeting 
point of the asymptotes at t ~ 1290 days lies at (87 ~ 0.244. 
At this point (87 for the fluid « 0.33, so the fluid is still 
moving at a significant fraction of the speed of light. 

According to eq. 1221 the predicted radius for the tran- 
sition to occur is i?NR ~ 0.38 parsec for the initial explosion 

For example, an Intel dual core 1600 MHz processor with 4 GB 
of ram. 
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Figure 1. ^7 at the shock front for the advanced EOS simulation 
(solid line), along with its asymptotic behaviour both at early and 
late stages. For com parison we have also plotted a prediction from 
iHuang et al] lll999h (see text). 

energy and circumburst density that we have used, corre- 
sponding to a lab frame time £nr ~ 450 days. We therefore 
conclude that the transition point from the relativistic to the 
nonrelativistic regime is far later than predicted by tivij. 

Also plo t ted in fig. [1] is the predicted value for f3j from 
IHuang et all (|l999T l. which we have implemented as follows. 
The starting point is 

d ^ - ^ ~ 1 (23) 
dm M e j + 27m ' 

the differential equation proposed by the authors to depict 
the expansion of GRB remnants, simplified to the adiabatic 
case. Here m denotes the rest mass of the swept-up medium 
and M e j the mass ejected from the GRB central engine. Our 
approach starting from the BM solution is a limiting case 
where M e j J. 0. The M e j term was included bv lHuang et al.l 
(1999) to incorporate a coasting phase. When solving eg. 1231 
we will use a very high 10 7 ) initial bulk fluid Lorentz 
factor 70 and by assuming M c j ^ E/2joc 2 we converge on 
the limiting scenario used in our simulations. Eq. [221 can be 
analytically solved to yield 

(7-l)Af cj c 2 + (7 2 -l)mc 2 = E, (24) 

which (numerically) leads to f(t) once we apply 

4 3 

m = gt-R womp, (25) 

and 

R(t) = f P{t)cAt. (26) 
Jo 

Here t is measured in the simulation lab frame (i.e. it does 
not refer to observer time). 

The resulting curve for (3y initially lies below the sim- 
ulation result, but ends up above at 4/3 times the simu- 
lation value. The initial and final slopes for the analytical 
(87 curve are c orrect by constructio n. We conclude that the 
approach from IHuang et all (|l999T i initially underestimates 
the BM phase and significantly overestimates the late stage 
flow velocity. The transition point between the relativistic 
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Figure 2. The resulting blast wave radii as a function of lab frame 
time for different simulations. The steady slope line shows the 
radius as predicted by the ST solution. The different simulations 
end up in the asymptotic regime with different radii: the T a d = 
5/3 ends up above the ST solution, the advanced EOS below the 
ST solution between the others and Y a ^ = 4/3 the lowest. The 
bottom curve shows the effective adiabatic index for the advanced 
EOS, minus 4/3. It starts at approximately zero at the left of the 
plot and proceeds to its asymptotic limit 1/3 in the nonrelativistic 
case. 



and nonrelativistic regime also lies at an earlier time for the 
analytical curve, closer to the analytically predicted 4nr- 



3.3.2 Blast wave radius 

In figure we plot the blast wave radius as a function of 
lab frame time for three different simulations: fixed adia- 
batic index at 4/3 and 5/3 and using the advanced EOS. 
Also we plot the radius as predicted from eq. l21l and. for the 
advanced EOS simulation, the difference between the effec- 
tive adiabatic index and its relativistic limit F a d = 4/3. The 
latter illustrates how relativistic the fluid still is in terms 
of temperature (as opposed to flow velocity). At the inter- 
secting point for the y/3 asymptotes at 1290 days the effec- 
tive adiabatic index is already quite close (r a d, e ff ~ 1.63) to 
its nonrelativistic limiting value. After 3800 days, when the 
time evolution of all the radii has become practically indis- 
tinguishable from R(t) oc t 2 ^ 5 from the ST solution we still 
see a difference between the different radii. At this time the 
ST radius is 1.358 parsec, the T a d = 4/3 radius is 1.197 par- 
sec, the T a d = 5/3 radius is 1.388 parsec and the advanced 
EOS radius 1.313 parsec. Taking the advanced EOS radius 
as a standard, this implies that ST overpredicts the radius 
at this stage by 3.4 percent, F a d = 5/3 overpredicts the ra- 
dius by 5.7 percent and T a d = 4/3 underpredicts the radius 
by 8.9 percent. Because all radii follow close to the same 
temporal evolution at this stage, these errors will only very 
gradually become smaller throughout the further evolution 
of the blast wave. Therefore, a derivation of the quantity 
E/po from the radius using the Sedov- Taylor equation l21l is 
likely to overpredict its value by approximately 18 percent 
within any time interval of practical interest. 
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Figure 3. Comoving number density profile. The profiles were 
taken at times corresponding to emission arrival times for the 
closest part of the shock front (i.e. with velocity directly towards 
the observer) at fO, 100, 1,000 and 10,000 days. Listed in increas- 
ing number and including the initial profile, these times corre- 
spond to lab frame times of 137, 387, 761, 2,227 and 12,583 days. 
Later times correspond to curves peaking further to the right in 
the plot. 
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Figure 4. Lab frame number density divided by 7 2 . This effec- 
tively scales the shock profile to 4 at the shock front. The same 
lab frame times as in fig. \3\ have been used. 



3.3.3 Density and energy profiles 

In figure [3] we have plotted the comoving fluid number den- 
sity profile (of the protons or the electrons, not both) at 
five different moments in time. Because later on we discuss 
spectra and light curves up to an observer time of 10,000 
days, we have chosen emission times corresponding to ar- 
rival times of the shock front (using t bs = t — R(t) /c) up to 
10,000 days as well. The earliest fluid profile shows the initial 
conditions calculated from the BM solution when the shock 
Lorentz factor is 10. After some time, the number density at 
the shock front can be seen to tend to the value predicted 
from the shock-jump conditions for a strong classical shock, 
which is no(r a( j + l)/(r ac j — 1) = 4no for the classical value 
of the adiabatic index 5/3. 

What is shown in figure [4] is an interesting feature of 
the blast wave, which is that the lab frame density directly 
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Figure 5. Thermal energy density profile for the same lab frame 
times as in fig. [3] 



behind the shock divided by the squared fluid Lorentz factor 
directly behind the shock, D/j 2 — 4po throughout the entire 
simulation. This can be seen analytically to hold from the 
Rankine-Hugoniot relations in both the ultrarelativistic and 
nonrelativistic case, even though the adiabatic indices have 
different fixed values, from 



d _ r ad + 1/ 7 



pot 2 r ad -i ' (27) 

which can be viewed as the relativistic generalization of the 
classical compression ratio and holds for arbitrary 7. When 
we use an advanced EOS, where we let r a d smoothly evolve 
from 4/3 to 5/3, we see from the figure that this generalized 
compression ratio remains very close to four even at inter- 
mediate times. We make use of this feature for the shock 
detection algorithm (see appendix I A2[) . 

In figure[S]we have plotted the thermal energy density at 
the same times as the number density. Unlike the density at 
the shock front, the thermal energy density is not expected 
to tend to a fixed value. The ST solution instead predicts a 
steep decline oc R(t)~ 3 , which is why the final shock front 
thermal energy density is many orders of magnitude smaller 
than the initial shock front thermal energy density. 



3.3.4 Magnetic field and particle acceleration 

We now turn to those quantities calculated in AMRVAC solely 
to aid in the construction of spectra and light curves, and 
that have no feedback on the dynamics. In figure [5] we have 
plotted 6b. Because we assumed the number of field lines 
through a fluid surface element to remain frozen (see eqn. 
I12p . the magnetic field energy density declined less rapidly 
than the thermal energy and as a consequence the local frac- 
tion eb increased. A discussion on the merit of our assump- 
tion about the magnetic field behaviour is outside the scope 
of this work (and from particle-in-cell simulatio ns it can cer- 
tainly be argued that it is not perfect, see e.g. IChang et al.l 
2008). However, our plot does show that at least it does not 
lead to unphysical values or strong inconsistencies. The max- 
imum value for eb found in fig.[6]is 0.037 (up from 0.01 at the 
shock front), which is not unreasonably large and, besides, 
occurs far downstream in a region that will contribute negli- 
gibly to the observed flux. We emphasize that ts is a relative 
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Figure 6. The fraction of the thermal energy that resides in the 
magnetic field energy eb f° r the same lab frame times as in figure 

El 



Figure 8. Upper cut-off Lorentz factor of the power law distri- 
bution 7j n , normalized to 1 at the shock front, for the same lab 
frame times as in fig. [3] 




Figure 7. £ni the fraction of electrons accelerated to a power law 
distribution for the same lab frame times as in fig. [3] 



measure and that both thermal and magnetic energy densi- 
ties drop steeply, both with respect to earlier times and with 
respect to their value at the shock front at any time. The 
numerical method presented in this paper for parametrizing 
the magnetic field energy density is quite general and can 
be readily modified to study different parametrizations. 

In figure [7] we have plotted the fraction £n of electrons 
that are accelerated to a power-law distribution. This frac- 
tion was taken to smoothly decrease from unity in the rela- 
tivistic regime down to 0.1 for our simulation in the nonrel- 
ativistic regime. The rightmost profile, with the shock front 
arriving at 10,000 days, has £n down to 0.16. 

In figure [8] we have plotted the normalized values for 
7m, the upper cut-off Lorentz factor of the power-law par- 
ticle distribution. Although formally 7 M should be reset to 
infinity at the shock front, we picked a value corresponding 
to a cut-off above 10 18 Hz through the entire simulation, for 
a fluid element heading directly towards the observer (see 
also section [231 and appendix^} . In our simulation settings, 
this results in 7 M peak values of the order 10 7 , and because 
these values are arbitrary as long as they are sufficiently 
large, we have normalized the 7 M profiles. The profiles show 



two things. First, they show the steep decline directly after 
the injection of new hot electrons. This steep decline and 
the ^ 7 orders or magnitude difference between shocked and 
unshocked 7 M are numerically challenging, which is why we 
implemented the logarithm of 7 M in our code instead (again 
see appendix |A"| . Second, the width of the profile is a mea- 
sure for the size of the hot region discussed in section 12.21 
It can be seen from the figure that the width of the pro- 
file increases over time. The width will nevertheless remain 
smaller than the width of the density profile by far. In our 
simulations, we resolve the 7 M profile and use it to determine 
the local refinement level. 



4 SPECTRA AND LIGHT CURVES 

Using the simulation data described in the previous section 
we have calculated spectra and light curves at various ob- 
servation times and frequencies. We have saved a total num- 
ber of 10,000 snapshots of the fluid profile, with 10.8 • 10 4 
seconds between consecutive snapshots, corresponding to a 
resolution cdt ^ 3 • 10 15 cm. Although this resolution is of 
the same order as the initial shock width, it is still sufficient 
at the early stage because the shock initially nearly keeps 
up with its own radiation. The effective resolution is given 
by cdt/T 2 - 3 • 10 13 cm, which is only a factor 10 larger 
than the spatial resolution of 10 12 cm and corresponds to 
a temporal resolution of At ^ 1000 seconds. It is therefore 
ensured that the blast wave in the initial stage is covered by 
over a hundred snapshots. 



4.1 Expected spectral and temporal behaviour 

The scaling behaviour for the critical frequencies and the 
flux is well known from analytical estimations assuming 
a homogeneous radiating slab directly behind the shock 
front with fluid proper ties determined via either the BM 
or ST solution (see e .g. Wiiers fc Galamalll999l ; iFrail et all 
l200d : iGranot fc Sari l2002h . We summarize the interstel- 
lar medium (ISM) scalings below, with denoting the 
peak frequency, va the synchrotron self-absorption critical 
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frequency and v c the cooling break frequency. At the ob- 
server times and frequencies in this paper we find either 
va < v m < v c or v m < va < Ve- 
in the relativistic limit, the corresponding scalings are 

v A < v m 

" A X 1 i -(3p+2)/2(p+4) ,,. ^ „ \*°) 



VA > V n 



V, 



oc t 

Vr OC t 



-3/2 
-1/2 



(29) 
(30) 



for the critical frequencies. Note that t now refers to observer 
times. The flux above both peak and self-absorption break 
scales as 

/ (l- P )/2 t 3(l- P )/4 i y<Vc 
,- P /2 t (2-3 P )/4 „>„„ " 



F oc 



(31) 



If va < Vm we get the following flux scaling below the peak 
break: 



F oc 



2,1/2 



uH 



1/3,1/2 



V < VA 

V > VA 



If va > Vm we have below the self-absorption break 



F oc 



1/ > v n 



In the nonrelativistic limit the scalings are 



t°'°, V A < V n 

VK X £-(3p-2)/(p+4) 



VA > V n 



Vm OC t 
V c OC t~ 



1/5 



(32) 



(33) 



(34) 

(35) 
(36) 



for the critical frequencies. The flux above both peak and 
self-absorption break scales as 



F oc 



Iy (l- J ,)/2 t (21-lS P )/10 ) l/<Vc 
u - P /2 t (^: ip)/ 2 v> 



(37) 



If va < Vm we get the following flux scaling below the peak 
break: 



F oc 



2,-2/5 



vH 



V < V A 

V > V A 



If va > Vm we have below the self-absorption break 

V 2 t 13/5 , V < Vm 



F oc 



i/ 5/2 f U/10 



V > V n 



(38) 



(39) 



The summary above shows that only the temporal be- 
haviour of the break frequencies and fluxes is altered by the 
transition to the nonrelativistic regime. We therefore do not 
expect spectra calculated from our simulations covering the 
transition to differ in slope from the slopes calculated above. 
The light curve slope, however, may differ. 

4.2 Spectra 

In figure [9] we have plotted spectra for a number of differ- 
ent observation times, ranging from 1 day to 10,000 days. 
For comparison we have also p lotted the different po wer law 
slopes for 1 day as predicted bv lGranot fc Saril ((2002) , where 
we have added a dependency on £n- We plot predictions for 
both £n = 1 and £n = 0.1. It can be seen that the simu- 
lated spectrum still lies closer to the £n = 1 prediction, just 
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Figure 9. Spectra at different observer times. The smooth curves 
show simulated spectra at different observer times: 1, 10, 100, 
1,000 and 10,000 days, with later observed spectra having lower 
flux in the high frequency range. For comparison we have included 
predicted slopes at the different power law regimes after 1 day, 
for both £n = 1 (solid line) and £n = 0.1 (dashed line). 



as we would expect for an early time spectrum. Because of 
shifts in both flux level and position of the spectral breaks 
for different values of £n, the flux does not always lie in be- 
tween the analytical predictions. For example, because the 
peak frequency Vm for the simulation lies close to that of 
the £n = 1 prediction and flux lies also closer to £n = 1 but 
below, the resulting flux at higher frequencies ends up below 
both predictions. 

Figure proves that our method works and that the 
asymptotic behaviour for the spectral slopes matches the 
predicted slopes. For frequencies above the self-absorption 
break and below the cooling break this merely confirms that 
the synchrotron spectral function Q(v' /v'm) has been imple- 
mented correctly. The flux at frequencies above the cooling 
break however, shows the consequence of a finite and evolv- 
ing upper cut-off 7m ■ A slope is reproduced that matches the 
prediction. It has been explained above and in EW09 how 
this slope now arises as a product of the interplay between 
the hot region and the blast wave width. 

At the low frequencies, where synchrotron self- 
absorption plays a role, the simulations also reproduce a 
spectral slope that corresponds to what was expected from 
analytical calculations. The flux level is now dictated by 
the radiative transfer equation through a medium that is no 
longer completely transparent at these frequencies. As dis- 
cussed in section |2"31 the resu l ting fl ux will differ by a factor 
of a few from iGranot fc Saril (|2002l ). due to a difference in 
approach when calculating the absorption coefficient from 
the particle distribution. 

We emphasize that fig [9] covers 10 orders of magnitude 
in frequency, 8 orders of magnitude in flux and four orders 
of magnitude in observer time. As we expected from analyt- 
ical calculations, the spectral slopes in the different power 
law regimes do not change over time. The transitions be- 
tween the different regimes are smooth. An explicit calcula- 
tion of the sharpnesses of the transitions will be presented 
in a follow-up paper. 
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Figure 10. Comparison of optical (at t ■ 10 14: Hz light curves for 
different equations of state. The top curve has T a( j = 4/3, the 
center curve the advanced EOS and the bottom curve r ad = 5/3. 
For clarity as few complications as possible are included: cooling 
and self-absorption are switched off, is fixed at 0.01 and £n = 
0.1 everywhere. Also plotted are the expected relativistic slope 
3(1 — p)/4 and nonrclativistic slope (21 — 15p)/10. 
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Figure 11. Direct comparison between thermal energy density 
e th P ron l cs f° r the different equations of state. The top profile has 
r ad = 4/3, the center curve the advanced EOS and the bottom 
profile has r ad = 5/3. All snapshots are taken at 515 days simu- 
lation time. The difference in radius between the blast waves is 2 
percent. 



4.3 Light curves 

We will use optical light curves to illustrate the consequences 
of the different assumptions and model parameters. In fig. 1101 
we present simulated light curves for simulations that differ 
only in the EOS used. Electron cooling and self-absorption 
have been disabled, es is fixed at 0.01 and £n = 0.1 ev- 
erywhere throughout the simulation. This allows for a clear 
view on both the effect of the EOS and of the transrelativis- 
tic break. The latter can be found at ^ 1000 days for all 
three simulations. This is somewhat earlier than the tran- 
sition time determined from the fluid flow in section 13.3.11 
which we determined to be around 1290 days (the difference 
is due to relativistic beaming). The transition time is also 
later than what is usually assumed for the nonrelativistic 
transition by nearly a factor three. 

The difference in flux from the different EOS assump- 
tions can be traced to the different thermal energy profiles 
(and hence, for fixed eb, to magnetic field energies that dif- 
fer with the same ratios), with r a( j =4/3 having the high- 
est e' th . This is illustrated in fig. 1111 The difference in peak 
thermal energy densities between the fixed adiabatic index 
simulations is a factor of two, as expected from the ratio 
(5/3 — l)/(4/3 — 1). Because the flux depends on the ther- 
mal energy via the magnetic field strength and 7 m (see equa- 
tions [5] and 0, the flux for r a( j = 4/3 is higher than that 
for T a d = 5/3. The light curve for the advanced equation of 
state lies between the two limiting cases, starting close to 
the 4/3 curve but moving to the 5/3 as the flow becomes 
nonrelativistic. This additional decrease in flux has the con- 
sequence that the advanced equation of state light curve will 
be slightly steeper in the transrelativistic phase than the fixed 
adiabatic index light curves. 

In figure [12] we show the effects of the detailed evolu- 
tion calculation of £n- Aside from the full simulations, we 
also perform two simulations that keep £n fixed throughout 
at either 1 or 0.1 but are otherwise identical to the full simu- 
lation. At early times in the radio, before the peak frequency 



0.05 -\„ 




1 10 100 1000 10000 

observer time (days) 



Figure 13. Fractional difference between complete and fixed 
£B = 0.01 simulation light curves. Solid line for radio, dashed 
line for optical. 



has passed, the £n = 1 curve lies above the full simulation 
curve, where at early times in the optical it lies below. 

Figure [13] shows the fractional difference between com- 
plete and fixed e_e = 0.01 simulation light curves (the 
light curves themselves lie very close to each other on a 
plot using a logarithmic scale), calculated via (infixed — 
Complete) /-^complete, where F is the flux. The figure shows 
that the late time light curves for the fixed e_e end up be- 
low the light curves that trace the evolution of the magnetic 
field. This can be understood from the fact that evolving 
e' B according to e' B oc (p') 4 ^ 3 implies a relative rise of the 
magnetic field energy density relative to (but still far below) 
the thermal energy when the flow becomes nonrelativistic. 
Fixing e_e forces the magnetic field energy density to follow 
p rad . The flux spans many orders of magnitude over time. 

A quantitative comparison between the slopes from the 
radio and optical light curves for the full simulations is 
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Figure 12. Left: Comparison between complete simulation light curve (solid line) at radio frequency 4.8 • 10 9 Hz and simulation curves 
where £n is kept fixed at 1.0 (dashed line) and 0.1 (dotted line) throughout. The complete curves start out close to £n = 1 but slowly 
evolve towards £n = 0.1. Right: same as left, only now for optical frequency 5 ■ 10 14 Hz. As with the radio light curves, the full curve 
starts near £n = 1 but turns to £n = 0.1. 
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Figure 14. Power law behaviour of the optical and radio light 
curves. The lines plot a, assuming for every two consecutive data 
points the relationship Fi+i = Fi(ti+i/ti) a . The solid line refers 
to the radio light curve and the dashed line to the optical light 
curve. The horizontal lines denote 0.5, 3(1 — p)/4 = —1.125, (2 — 
3p)/4 = -1.375, ((21-15p)/10) = -1.65, (4-3p)/2 = 1.75 from 
top to bottom. 



shown in figure 1141 The horizontal lines in the plot indi- 
cate expected asymptotic values for the power law scalings. 
In the relativistic limit, the expected slope is 1/2 before pas- 
sage of z/ m and —1.125 after (using p — 2.5). After the cool- 
ing break passes, a further steepening to —1.375 is expected. 
In the nonrelativistic regime the expected slopes before and 
after passage of the cooling break are —1.65 and —1.75 re- 
spectively. The plot shows that the relativistic slopes are 
matched very well. The radio light curve quickly tends to 
1/2 and after passage of the peak frequency it moves in 
^ 95 days to —1.125, where it remains until the onset of the 
nonrelativistic break time. The optical light curve starts out 
in the intermediate regime from the passage of v m , with the 
passage of the cooling break coming too early for the light 
curve to settle into the pre-cooling break slope of —1.125. 



The post-cooling break slope —1.375 is obtained instead and 
is again maintained until the nonrelativistic break. The light 
curve slopes in the nonrelativistic regime are less steep than 
expected. A number of factors play a role here, as discussed 
above. The advanced EOS leads to a steepening of the de- 
cay during the transition phase (which lasts well over 10,000 
days), whereas the increase in e' B relative to e' th and the de- 
crease in £n (leading to an increase in energy per particle) 
lead to less steep decay. The change in slope in the nonrel- 
ativistic regime is the result of the interplay between these 
different factors, with the end result being a slope less steep 
than expected. The final nonrelativistic slopes differ signif- 
icantly from those expected from analytical models, and this 
has a large impact on fitting models to observational data. 



5 GRB030329 

In the preceding section we have systematically explored the 
different aspects of transrelativistic blast wave afterglows 
with respect to dynamics and radiation for standard values 
of the input parameters. We now qualitatively compare ra- 
dio data for GRB030329 to simulation results using physical 
parameters for this GRB established by earlier authors as 
input. GRB030329 is one of the closest and brightest GRBs 
for which an afterglow was found. Because of this brightness, 
the afterglow could be monitored for an extended period of 
time at various wavelengt hs and after six year s its radio sig- 
nal is still being observed (|Kamble et al.ll2009h . GRB030329 
is a good example to use to illustrate the various aspects of 
the radiation code. 

The redshift of GRB030329 has been determined to 
be 2 = 0.1685 l|Greiner et al.l 120031 ). which leads to a 
luminosity distance of 2.4747 ■ 10 27 cm (for a flat uni- 
verse with Q M = 0.27, Qa = 0.73 and H = 71 km 
s _1 Mpc" 1 ). Various authors have determined the phys- 
ical properties of the GRB from ana l ytical model f its to 
the data (e g. IWillingale et all 120041: " 



Berger et all 120031: 



Sheth et all 120031: IVan der Horst et all 120051 : iHuang et all 
20061 : Van der Horst et al. 20081 ) with various assumptions 
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for the jet str ucture. Here we take t he ph ysical parameters 
established bv lVan der Horst et al.l l|2008l ). From their con- 
clusion for the jet break time and cooling frequency at this 
time, and assuming equipartition between accelerated parti- 
cle energy and magnetic field energy (i.e. a fixed = eg in 
their model), we arrive at E = 2.6 • 10 51 erg (for a spherical 
explosion), n = 0.78 cm -3 , p = 2.1, e E — e B — 0.27. We 
assume a homogeneous medium and we set the hydrogen 
mass fraction in this medium to unity. IVan der Horst et al.l 
(2008) fix £jv at unity, but we use a nonrelativistic limit 
Civ.NR = 0.1. Because GRB030329 shows clear evidence of 
a collimated outflow, it is no longer sufficient to assume a 
spherical explosion. When calculating emission from a jet, 
we assume a hard-edged jet with opening angle 22 degrees 
and no lateral spreading. 

We have plotted light curves at 15 GHz, 4.8 GHz and 
1.4 GHz in figure [151 including data points from the West- 
erbork S ynthesis Radio Telescope (WSRT) (4.8 GHz and 
1.4 GHz. lVan der Horst et al. 2005) for comparison a nd the 
Very Large Array (VLA) (15 GHz. feerger et alj|200j ). Two 
things are clearly visible. First, our simulated light curves 
still differ strongly from the data, although largely the same 
input parameters have been used for the blast wave simula- 
tions as those that were derived from fitting to the dataset 
using an analy tical model for the blast wa ve. The different 
assumptions in IVan der Horst et all (|2005h account for this 
in part, but nevertheless this demonstrates once more the 
need for detailed fit prescriptions from simulations (a simi- 
lar conclusion was drawn in EW09 for the ultra-relativistic 
case). Second, the counterjet contribution will stand out 
clearly for a hard-edged jet model. For now, the compari- 
son between simulation and data is still qualitative. Newer 
data are available and once the simulation input parameters 
are fine-tuned with respect to the data as well (as opposed 
to estimated using an analytical fit to the data), it should 
be possible to address the rise of the counterjet in a more 
quantitative fashion. 

Because of the equipartion constraint on tj and tE, 
both were given a relatively high value of 0.27 at the shock 
front. In the nonrelativistic regime, the magnetic energy den- 
sity will grow relative to the thermal energy density further 
downstream (although both will decrease strongly in abso- 
lute value). At t e ^ 34.7 yrs, the last time covered by our 
simulation (set up to cover 10,000 days in observer time) 
we find that tB has risen to approximately 0.36 at the back 
of the blast wave (0.43 where the blast wave density again 
equals the upstream density). Even further downstream, 
when the density has fallen three orders of magnitude below 
the upstream density, es peaks at 1.28. This is not unphys- 
ical, but merely an indication that magnetic fields have be- 
come dynamically important in a region of the fluid which 
has no consequence for the light curve. 

5.1 Low frequency array light curves and resolved 
images 

Figure [16] shows predicted light curves at the very low fre- 
quencies that can be explored in the near future by ra- 
dio telescopes such as the Low Frequency Array (LOFAR), 
assuming four hours o f integration time, 25 core station s 
and 25 remote stations (jNiiboer fc Pandev- Pommiersll2009h . 
GRBs are among the prime targets for LOFAR's Transient 



Key Project <|Fenderll2006h . Most of the time all light curves 
lie below the self-absorption break. This, in combination 
with the Vm break, a hard-edged jet model and the turnover 
to the nonrelativistic regime, leads to an interesting double 
peak structure of the light curve. First the signal rises, ac- 
cording to the relativistic rise in the self-absorption regime 
that predicts a slope of 1/2. After ^ 4 days a clear jet break 
is seen and the resulting drop in slope leads to a decreasing 
signal again. Around circa 150 days the critical frequency 
Vm passes through the observed frequency band. The slope 
of the spherical explosion changes accordingly towards the 
predicted relativistic 5/4. Around approximately 600 days 
the blast wave has become nonrelativistic and the counterjet 
starts to contribute (but is still overwhelmed by the forward 
jet). The predicted nonrelativistic slope for the spherical ex- 
plosion is now 11/10. 

We have included LOFAR detection thresholds for four 
hours of integration ti me. These sensitivity l i mits a re higher 
than those presented in IVan der Horst et all l|2008h , because 
LOFAR has been scaled down in the meantime. The spher- 
ical explosion energy is an overestimation of the actual ex- 
plosion energy and the flux levels corresponding to the jet 
simulations lie closer to what will actually be received. How- 
ever, from fig. [15] it is clear that our qualitative compar- 
ison systematically underestimates the actual flux levels. 
Also, the integration time used in LOFAR can easily be 
increased, even up to days. Fig. 1 161 therefore does not mean 
that GRB030329 will not be observable by LOFAR, but only 
that a larger integration time than four hours is likely re- 
quired. 

For 200 MHz we have calculated spatially resolved im- 
ages as well, for spherical explosions. Three images are pre- 
sented in figure [T7l for three different observer times. They 
show three qualitatively different types of behaviour. At 15 
days a limb-brightened image is observed, whereas at 240 
days the image on the sky becomes limb-darkened. At 3900 
days another structure is visible and a brighter ring exists 
within the image, at a radius of ^ 10 18 cm. This is a result 
of the self absorption break va being different for different 
emitting regions of the blast w ave. Th e se im ages are fully 
consistent with predictions from lGranotl (|2007f l for the ultra- 
relativistic case. 



6 SUMMARY AND CONCLUSIONS 

In this paper we present the results of detailed dynami- 
cal simulations of GRB afterglow blast waves decelerating 
from relativistic to nonrelativistic speeds, as well as spec- 
tra and light curves calculate d from these simulations us - 
ing a method first described in I Van Eerten fc Wiiersl (|2009h 
(EW09) that we have extended to include more details of 
synchrotron radiation. We summarize our results and con- 
clusions below. 

We have performed, for the first time, hydrodynami- 
cal simulations of decelerating relativistic blast waves using 
adaptive mesh refinement techniques including a parametri- 
sation for a shock accelerated electron distribution radiating 
via synchrotron radiation. From these simulated blast waves 
we have calculated light curves and spectra at various ob- 
server times and frequencies. An advanced equation of state 
was used for the dynamical simulations, with an effective 
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Figure 15. Left: Light curves at 15 Ghz, plus data. The upper curve is calculated from a spherical explosion, the bottom from a hard- 
edged jet with opening angle 22 degrees, while the crosses are the data. The receding jet is clearly visible. The bottom right curve is the 
radiation from the counterjet alone. The transition to the nonrelativistic regime can be seen from the spherical explosion simulation at 
around the same time when the counterjet becomes visible. Center: Light curves at 4.8 Ghz, plus data. Right: Light curves at 1.4 Ghz, 
plus data. 
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Figure 16. Left: Simulated light curve at 200 MHz for GRB030329, top curve for spherical explosion and bottom curve for hard-edged 
jet with opening angle 22 degrees. We have drawn the following slopes from left to right: 1/2, 5/4 and 11/10. LOFAR sensitivity for 25 
core and 25 remote stations after four hours of integration time is 0.273 mjy and indicated by the horizontal line. Center: Simulated light 
curve at 120 MHz for GRB030329. LOFAR sensitivity is 0.145 mjy. Right: Simulated light curve at 75 MHz for GRB030329. LOFAR 
sensitivity is 4.2 mjy, too high to be shown in the plot. 




Figure 17. Radio images at 200 MHz. Left: after ^ 16 days. The intensity increases monotonically outward. The outer radius is 1.91 10 17 
cm. Center: after ^ 270 days. The intensity decreases monotonically outward. The outer radius is 9.7- 10 17 cm. Right: after ^ 4300 days. 
A central bright ring with radius ^ 10 18 cm appears. At larger radii the intensity decreases monotonically. The outer radius is 3.4 ■ 10 18 
cm. 
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adiabatic index smoothly varying between the relativistic 
and nonrelativistic limit. Three additional parameters were 
traced during hydrodynamical evolution: maximum acceler- 
ated particle Lorentz factor, magnetic field energy density 
and accelerated particles number density. We assumed that 
fewer particles were accelerated by shocks that are less rel- 
ativistic. To obtain the observed flux including synchrotron 
self-absorption, a set of linear transfer equations were solved 
for beams traversing through the blast wave. This method 
expands upon EW09 by including self-absorption and dy- 
namically calculated electron cooling. 

We have used standard assumptions for the GRB explo- 
sion energy 10 52 erg) and circumburst particle number 
density 1 cm -3 ) for a homogeneous medium and parti- 
cle acceleration and magnetic field parameters. By directly 
comparing against various analytical models and expected 
limiting behaviour, we draw a number of conclusions about 
the dynamics of our simulations: 

• We find that the transition of /3y directly behind the 
shock front from the relativistic to the nonrelativistic regime 
occurs later than expected, around ^ 1290 days rather than 

450 days, for the standard model parameters. 

• An an al ytical calculation of /3y according to 
iHuang et all jl999h is found to overestimate the late 
time values by a factor 4/3. 

• Directly applying the Sedov- Taylor solution to late time 
afterglow evolution is found to overestimate the radius by a 
few percent and keeping the adiabatic index fixed through- 
out the evolution of the blast wave will lead to systematic 
differences of as much as ten percent. 

• The density jump across the shock may be arbitrarily 
high for relativistic shocks, but will be a factor of four in 
the nonrelativistic regime. This is known from the shock 
jump conditions. Our simulations show that the quantity 
D/j 2 , a combination of lab frame density and Lorentz factor 
directly behind the shock, will remain close to four times 
the unshocked density throughout the entire simulation, even 
though the effective adiabatic index evolves from relativistic 
to nonrelativistic. 

• If we assume the number of magnetic field lines through 
the surface of a fluid element a constant, the magnetic field 
energy will become relatively larger compared to the thermal 
energy. It will remain a small fraction however (assuming 
only a small amount of energy is used for magnetic field cre- 
ation accross the shock). Our approach allows for different 
assumptions on the magnetic field energy evolution. 

• The upper cut-off Lorentz factor 7^ for the shock- 
accelerated relativistic power law electron distribution de- 
creases on a distance scale much smaller than the width of 
the blast wave due to synchrotron losses and determines the 
shape of the spectrum near and above the cooling break. 

Using the output from the dynamical simulations, we 
calculate the flux. The following general conclusions are 
drawn for the radiation: 

• Calculated light curves show a transition between the 
relativistic and nonrelativistic regime at around 1000 days 
in observer time, again later than expected. 

• The observed fluxes for different assumptions on the 
equation of state may differ by a factor of a few. This is a 
direct consequence of the amount of thermal energy (and 



therefore magnetic field energy) directly behind the shock 
front. 

• Implementing a changing effective adiabatic index has 
the consequence that the resulting light curve will slowly 
evolve from the relativistic limiting value to the nonrela- 
tivistic value. This transition takes tens of years in observer 
time and will lead to a steeper decay in the afterglow light 
curve than predicted by analytical models assuming a fixed 
index. 

• This steepening is a smaller effect than the combined 
effect of evolving the magnetic energy density and the ac- 
celerated particle number density. When all effects are in- 
cluded, the final light curve slopes differ markedly from the 
analytically expected values. This implies a significant com- 
plication for late time afterglow modeling. 

We have applied our approach to GRB030 329 as well, 
using physics parameters derived by IVan der Horst et al.l 
(2008) using an analytical model. It is shown that the re- 
sulting radio light curves differ up to an order of magnitude 
between simulation and analytical model, although this can 
be partly attributed to some different assumptions. Assum- 
ing a hard edged jet with an opening angle of 22 degrees, 
our simulated light curves show a rebrightening due to the 
counterjet around 1000 days. Simulated curves at radio fre- 
quencies that will be observable using LOFAR show that 
four hours of integration time is likely not sufficient to dis- 
tinguish the signal from the noise and a larger integration 
time is required. Finally, spatially resolved images show a 
bright ring that, depending on the precise power law regime 
that is observed, may be located not only in the center or on 
the edge but also at intermediate radii withi n the afterglow 
image. This is consistent with earlier work bv lGranotl (|2007T l 
on afterglow images i n the relativistic phase . 

A recent paper (|Zhang fc MacFadvenl l2009h has ap- 
peared discussing afterglow blast waves decelerating to non- 
relativistic velocities using twodimensional simulations. The 
authors find that lateral expansion of a relativistic GRB 
jet is a very slow process and that the jet break is mostly 
due to the edges of the jet becoming visible. This im- 
plies the hard edged jet model that we have applied to 
GRB030329 is sufficient to model the jet break at ^ 4 days. 
IZhang fc MacFadvenl (|200gh do not include synchrotron self- 
absorption and calculate the cooling break by assuming the 
cooling time throughout the entire blast wave equal to the 
grid time. 

The approach to calculating light curves and spectra 
from generic fluid simulations that we present in this paper 
assumes that synchrotron radiation is the dominant radia- 
tive process, that particle acceleration takes place in a region 
far smaller than the blast wave width and that the feedback 
on the dynamics from the radiation is negligible. We briefly 
address these issues in appendix [Bl 
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APPENDIX A: NUMERICAL 
IMPLEMENTATION 

Al Partial differential equations 

AMRVAC was written to solve a system of coupled partial 
differential equations. When adding additional equations to 
the solver, it is therefore best to use partial differential equa- 
tions. In the case of the magnetic field energy e' B we start 
by rewriting equation 1131 as 



d 



d 



dt P 'i/i + v, d^^-°- (A1) 

If we multiply this equation by p'7 and add to this the con- 
tinuity equation 



d 



-jpjv =0, 
which we first multiply by e' B /p' 4 ^ 3 , we obtain 



d 'ye' B d 'ye'sv* 



+ 



n'l/3 



= 0. 



(A2) 



(A3) 



dt p'Vs 1 dx l p' 

This is the type of conservation equation that ARMVAC is 
specialized in, and it is therefore the quantity ^J 3 that we 
calculate in AMRVAC. 

For the evolution of the upper cut-off 7 M we follow a 
similar procedure. We start by simplifying equation [9] to 



d p 



'1/3 



dt 7m 



1 

7 



(A4) 



where a = or /Qnm c c, and t refers to lab frame time (i.e. 
emission time). The Lorentz factor in the source term arises 
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Figure Al. Normalized profiles of 7^ (solid line and squares), p' 
(dashed line and circles) and §N (dotted line and triangles). The 
actual values at the shock front arc 1.9 • 10 7 , 7.0 ■ 1CT 24 g cm J 
and 0.31 respectively. The simulation time is 1400 days. In the 
lower left corner we zoom in on the shock front, showing exactly 
where we reset 7^ and £n. 



the reset values of n^ cc and e' B depend on the fluid variables 
directly behind the shock front and it is therefore important 
that we determine the position of the shock front as accu- 
rately as possible. Mathematically speaking, a shock is a 
discontinuity in the flow variables with a sudden increase in 
entropy across the discontinuity. In practice, however, find- 
ing a shock in a numerical approximation is more involved, 
both due to numerical shock diffusion and because, strictly 
speaking, there is a shock discontinuity across every grid cell 
boundary. 

This has the consequence that if we try to find shocks 
by checking for discontinuities or for entropy jumps, we will 
find both shocks all over the numerical diffused shock region 
and at a random variety of positions where the numerical 
noise happens to rise above a predetermined shock treshold. 
This then implies that we keep on resetting the additional 
quantities over some region, something which is especially 
unwanted in the case of 7 M , given our approach where we 
take a fluid cell to contain a collection of electrons that have 
been shocked exactly at the same time and we critically rely 
on the size of the hot region (see section 12.21 and EW09, 
appendix D, for details). 



when we write the comoving time derivative in the lab frame. 
We can now follow a procedure similar to what we did for 
the magnetic energy density, but first we rewrite the equa- 
tion above once more for numerical reasons. The quantity 
7m varies over many orders of magnitude in a very short 
time span and any quantity that depends on 7 M linearly is 
therefore difficult to deal with numerically. We solve this by 
rewriting equation I A4I into 

i„C^. (A5) 
dt 7m 7 V ; 

Although there still is a linear dependence on 7 M in the 
source term, in practice this equation provides a better start- 
ing point for AMRVAC. From combining with the continuity 
equation we get 

d ,. p' 1/3 d 4 , p' 1 / 3 , , R , 2 .... 
_ 7p ln __ + _ v lp m _ = aw B , (A6) 

with 7p' ln p' 1//3 /7m the quantity of interest. A similar ap- 
proach to tracing the effect of cooling in the context of rel- 
ativist ic blast waves has also been taken by iDownes et al.l 
121)1)2). Although in our formalism 7 M at the shock front 
should be reset to infinity, and therefore -yp' bxp n ' 3 /j^ to 
minus infinity, we just take a very low value for I/7J1/ in 
order to minimize numerical diffusion. In our simulations, 
this arbitrarily low value corresponds to a hard cut-off of 
the spectrum above v ^ 10 ls Hz, at frequencies sufficiently 
far above our observation range to be of no consequence. 
The 'real' 7 M catches up with the numerical 7^ almost in- 
stantaneously. 

A2 Shock detection method 

In total, AMRVAC now calculates the evolution of three addi- 
tional quantities: n acc (using equation [TS]) , r ye' B /p' 1 '' 3 (using 
equation I A3[) and 7/9' ln(p' 1/ ' 3 /7 M ) (using equation I A6|) . All 
three quantities get reset wherever a shock is detected. Both 



Because the shocked particle number density and the 
magnetic field density depend directly on the fluid variables, 
using, for example, a j ump in /?7 as a trigger, as has been 
done bv lDownes et al.l(|2002l ). is not an option in these cases 
either. Although it serves as an excellent indicator of the 
front of a shock, it will not point us to a location where we 
can find information on the strength of the shock, but to an 
arbitrarily defined position just in front of that. 

In this paper we solve the issues of shock detection with 
two shock detection algorithms, both of them making use of 
the fact that -D/7 2 directly behind the shock is four times 
the density just in front of the shock. For n'^ cc and e' B we 
define the shock front to be at the peak of the Lorentz factor 
profile, in the region where D/j 2 > 3.5po- The numerical 
constant is arbitrary and could be taken closer to 4. With 
this method we ensure that the shock is detected at those 
positions where the fluid quantities are sufficiently close to 
their peak values, although multiple shock peaks may be 
detected in close proximity of each other due to numerical 
noise. 

For 7^ it is essential that we only detect a single shock 
front. Here we care less about the precise fluid variable val- 
ues. For the purpose of resetting 7^ we define the shock front 
to be at that position where D/poj 2 crosses the value 3.5. 
For a single shock front, this only happens once. Although, 
in principle, 7/3 could have been used instead of D/po7 2 , 
the latter offers the significant advantage that it does not 
change in scale over the course of the simulation and always 
remains close to 4, whereas 7/3 becomes arbitrarily small. 

Figure lATI illustrates the use of the Lorentz factor pro- 
file peak as a shock detector. It shows that the numerical 
diffusion is really very small and that 7^ changes over a 
significantly smaller spatial scale than p'. 
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A3 Synchrotron self-absorption 

Equation 1 181 can also be expressed as 



1m ; 

/v>' 2 / Ht>(£-) x 



(i - -^f - 4 

7m 



where 



K = C 



V3g*B' 
87rm?c 2 



(A7) 



(A8) 



Here we have used the fact that A e (7g) oc y' e (l — ^Hm) 
(i.e. a slightly modified powerlaw distribution). The scal- 
ing factor (C) of A/e(7e) is determined in terms of y' M and 
7m from the requirement that the total number of acceler- 
ated electrons constitutes a fixed fraction £jv of the available 
electrons. The symbol V denotes the pitch angle (the angle 
between magnetic field and particle velocity) averaged ver- 
sion of the synchrotron function, a dimensionless function 
representing the shape of the synchrotron spectrum for a 
single electron in the same way Q represents the spectrum 
of a distribution of particles. The v' cr ^ in the argument is 
connected to 7^ via equation [7] (see also EW09). 



By changing variables from j' e to y — 



we obtain: 



K 1 4ivm e c\-% ,-E±i 

v 2 X 



eCV 

v ) 



, K ft. 

Q "' 2 \3q e B> 

[h(VM,Vm,) + h(yM,Vm)], 

where the quantities Ii{yM,Hm) and l2(yM,ym) are: 



Vrm, 

h = ( P + 2) J dyP(y) 



y 



(A9) 



(A10) 



and 



h = (p-2)yF I dyr(y)y E ^(l~(^ 2 y- 3 . (All) 



As in the case of Q(yM,y m ) (see lVan Eerten fc Wiiersl 

2009) values of Ii(y M ,y m ) and Ii{yM,Vm) are tabulated 
for moderate yu and y m , whereas their limiting behavior, 
for extreme values of yu and y m is analytically estimated. 
Namely, if yulym 1 the integrals of both 7i and I2 reduce 
to the expression inside the integral, evaluated at y m , mul- 
tiplied by the appropriate range in y-space, i.e. (y m — yin)- 

For yu «C 1 the integrals' behaviour becomes hard to 
analytically estimate, especially for general values of y m and 
p. Instead, we fit approximate expressions to the values ex- 
trapolated from the tables. 

In the case that 2> 1 and y m is outside the tabulated 
values, we can break the integral into two parts by using the 
last tabulated value y m - For I2 the formula is: 
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Figure A2. Intensity and refinement levels perpendicular to the 
axis between the observer and the source. The maximum refine- 
ment level drops quickly to zero away from the edge of the jet. 
The intensity has been rescaled to an arbitrary scale suitable for 
direct comparison between intensity profile and refinement levels. 
Note that the lowest refinement level is zero. 
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For Ji the result is identical, only the terms inside the 
square brackets (including Q(x)) have to be evaluated for 
p -> p+ 1. 

Finally, for yu 3> 1 the result of both integrals is ap- 
proximated by zero. 



A4 Adaptive mesh and linear radiative transfer 

We do not integrate over A in equation [5] directly, but re- 
solve the different rays instead. After the integral over t e 
is finished (i.e. the bundle of linear radiative transfer equa- 
tions is solved) we can integrate over A to obtain the flux, 
while the unintegrated result provides a resolved picture of 
the emission from the fluid. For spherically symmetric fluid 
flow or for an observer positioned along the symmetry axis 
of a jet, the intensity on surface A is symmetric around a 
central point. Because the fluid itself moves at nearly the 
speed of light, it is not a priori clear how many rays need 
to be included and how they should be spaced along A in 
order to obtain a good resolution. An efficient response to 
this dilemma is to apply the adaptive-mesh refinement con- 
cept to A. The equidistant surface (EDS) A contains a grid 
with every grid cell containing the intermediate results for a 
single ray. Every four neighbouring cells in each direction on 
the EDS are grouped together in a single block. The EDS 
area dA that each cell represents may differ, and if the res- 
olution threatens to become too low to adequately capture 
the radiation profile, a block will be split in half along each 
direction, spawning new blocks that represent half the size 
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of the parent block along each direction. The refinement cri- 
terion that is used is that the combined flux from a given 
block must not differ by more than 1 percent (or a lower 
treshold, as set by the user) from the combined flux from a 
coarsened version of the block where only the odd cells are 
taken into account (with the odd cells representing an ap- 
propiately increased surface element). Neighbouring blocks 
may differ one refinement level at most. We have plotted 
an example of this strategy in figure IA2I In practice we set 
the maximum refinement level similar to that of the fluid 
simulation. We also use the fluid simulation grid refinement 
structure to determine the starting refinement structure of 
the EDS at each iteration for the transfer equation solver, in 
order to make sure that we will also capture the blast wave 
when it still has a small radius. 



APPENDIX B: 
MODEL 



APPLICABILITY OF OUR 



The radiation code is written to be generally applicable to 
output from relativistic fluid dynamics simulations. How- 
ever, a number of assumptions and simplifications have been 
made that are dependent on the physical context. In this ap- 
pendix we briefly discuss the consequences and relevance of 
our assumptions in the case of GRB afterglow blast waves 
decelerating down to nonrelativistic speeds. We discuss the 
relevance of an alternative radiative process, inverse Comp- 
ton scattering, of our assumption that particle acceleration 
takes places in a region much smaller than the blast wave 
width and of adiabatic expansion of the blast wave with the 
radiation losses having no effect on the dynamics. 



Bl Importance of inverse Compton scattering 

A limitation to the applicability of our approach arises from 
the fact that inverse Compton (IC) radiation, which is not 
calculated, becomes important when the ratio Ps yn /P{c a P~ 
proaches, or drops below unity. This ratio is also equal to the 
ratio between th e corresponding energy field s that power the 
emission e' B /e' ph <|Rvbicki fc Lightmajjl986h . with e' ph being 
the energy density of the (synchrotron) radiation field. The 
effect of IC emission o n the emitted spec tra has been thor- 
oughly investigated in lSari &P Esin (200l|). In this paper we 
focus only on its influence on cooling rates, as the high- 
energy synchrotron spectrum is expected to dominate IC 
emission for a wide range of physical parameters and radii. 

Instead, however, of calculating the entire photon en- 
ergy density due to synchrotron radiation we can use the 
fact that the cross-secti on for IC scattering drops fast be- 
yond the Thomson limit (|Blumenthal fc Gouldlll97Ch . Thus, 
we can define an 'effective' photon field for an electron of 
Lorentz factor 7^ as 
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where M^ hom 



m °, c h (with h denoting Planck's constant) 
is the photon frequency for which the scattering occurs 
marginally within the Thomson regime for a head-on col- 



lision and I' u i syn is the synchrotron specific intensity. Ap- 
proximating the specific intensity by 
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and employing the analytical relations of the BM solution 
we find that right behind the shock front 

§P » 7.5-10 16 f{p) (5- W-^e^e^ x 



4 U)— R^r*, (B3) 

where f(p) = (p - l) p ~ 2 (p - 2) 1 ~ p (3 - p), T is the Lorentz 
factor of the shock front and R the shock radius. By plugging 
in standard values of this paper (£n = 1, p = 2.5, no = 1 



10" 



€ E = 10 , E = 10 J erg) and making 



further use of the BM equations we find for 7^ 
1.8- 10" 10 i? - 5 . 



P' 
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This means that IC will dominate synchrotron energy losses 
for the lowest energy electrons throughout the relativis- 
tic phase of the fluid. A comparison of IC to adiabatic 
cooling, using the synchrotron loss term from eq. [9] and 
(d7m/dt) syn /(d7^/ dt)ic = P^/Pic, gives 
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(B5) 



(df m /dt) t0 

The corresponding radius after which adiabatic expansion 
will sharply take over is about 1.6 • 10 17 cm. Moreover, for 
an electron of energy 7c = 10 4 7^ (i.e. on the order of 7^), 
synchrotron losses will prevail at approximately 3 ■ 10 17 cm. 
Therefore it is only early on, and certainly not close to the 
subrelativistic transition, that IC cooling will affect the evo- 
lution of 7^ or 7m for the assumptions made in this paper. 

B2 Gyral radius 

The gyroradius for an electron with Lorentz factor 7 C is given 
by r'„ — "f c m e c 2 j q c B' . Using the BM solution and eq. [7] we 
find that for the most energetic electrons at the shock front 
that contribute to received flux within the frequency range 
under consideration (10 s — 10 18 Hz) this radius lies below 



= 5.6- 10 



.75 1/2 
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where i^ut-off is the cut-off frequency in the lab frame, used 
to set 7m at the shock front (i.e. a frequency safely above 
10 18 Hz). This should be compared against the size of the 
hot region, a measure of the spatial distance over which 
electrons cool significantly. The cooling time for high energy 
electrons that cool on a scale much shorter than the scale 
over which the fluid variables change, is approximately equal 
to t coo i w 67rm e c/<7T(-B') 2 7M when the electron cools down 
to 7 M . Using the self-similar parameter as an intermediate 
step, this can be linked to a spatial size of the hot region for 
7m: 



X o m-52 -1/2 D 33/8 

S hot =8-10 Kut_oS R cm - 



(B7) 

Thus, for ^ cu t-off = 10 18-21 Hz the gyroradius of the most en- 
ergetic electrons is many orders of magnitude smaller than 
the size of the corresponding hot region, justifying our as- 
sumption that particle acceleration takes place locally near 
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the shock front (which we have implemented by a local injec- 
tion of hot electrons) and the use of an advection equation 
to model the evolution of jm ■ 

B3 Feedback on the dynamics 

A last issue is the possibility of the radiative energy losses 
becoming comparable to the initial energy load of the fireball 
(E) that would imply a considerable impact on the dynamics 
of the flow. This could be explicitly quantified by calculat- 
ing the total radiative output during the simulations and 
comparing it to the explosion energy. However, we can ad- 
dress this issue in a more qualitative manner by noting that 
the low energy electrons cool predominantly by causing the 
expansion of the volume they are occupying (slow cooling), 
even at the shock front. Moreover, for values of p > 2 (which 
is the case under consideration) these electrons are the main 
energy carriers. In combination with the fact that the total 
energy residing in relativistic electrons is limited by £e (typ- 
ically of the order of 10%), we are confident that the total 
energy radiated through synchrotron, especially in the sub- 
relativistic regime, will be orders of magnitude smaller than 
E, and thus not affect considerably the dynamics. 



